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ZIP3D 

An Elastic-Plastic Finite-Element 
Program for Cracked Bodies 


I. SUMMARY 

ZIP3D is an elastic and an elastic-plastic finite-element program to 
analyze cracks in three-dimensional solids. The program may also be used to 
analyze uncracked bodies. For crack problems, the program has several 
unique features including the calculation of mixed-mode strain energy 
release rates using the three-dimensional virtual crack closure technique, 
the calculation of the J-integral using the equivalent domain integral 
method, the capability to extend the crack front under monotonic or cyclic 
loading, and the capability to close or open the crack surfaces during 
cyclic loading. This report includes three sections: a theoretical section, 
a user manual section, and two example problems (with input and output 
files). The theories behind the various aspects of the program are 
explained briefly in the theoretical section. Line-by-line data preparation 
is presented in the user manual. Input data and results for an elastic 
analysis of a surface crack in a plate and for an elastic-plastic analysis 
of a single-edge-crack-tension specimen are presented in the example 
section. 

II. THEORETICAL MANUAL 

INTRODUCTION 

ZIP3D is an advanced finite-element code developed to analyze cracks in 
elastic or elastic-plastic bodies. The stress-intensity factor (or strain- 
energy-release rate) for elastic materials and the J-integral for elastic- 
plastic materials are important parameters for damage tolerant design of 
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structures. For simple configurations and loading conditions, closed-form 
analytical techniques can be used to calculate the fracture mechanics 
parameters. However, for general structural configurations and loading 
conditions, and non-linear material behavior, the finite-element (FE) method 
is the best alternative. A large number of FE codes for both linear and 
nonlinear stress analyses of two- and three-dimensional bodies are available 
in the public domain and also through private companies. These codes are 
suitable for general engineering applications to calculate the usual 
quantities like stress, strain, displacements, forces, buckling loads or 
post buckling responses. Very few codes have the capability to calculate 
fracture mechanics parameters and those that are available are restricted to 
2-D (plane) problems or for evaluating only total strain-energy-release 
rates for three-dimensional (3-D) problems. 

The purpose of this report is to document a 3-D elastic-plastic finite- 
element program developed at the NASA Langley Research Center that 
calculates mixed-mode fracture mechanics parameters in cracked solids. 
However, bodies without cracks may also be analyzed to obtain stresses, 
strains and displacement fields. In this code, eight-node isoparametric 
elements and small -strain deformation theory are used. A special reduced- 
shear integration scheme [1] is provided for bending dominant problems. The 
elastic-plastic analysis is based on the incremental plasticity theory using 
the von Mises yield criterion, isotropic hardening, and Drucker's flow rule. 
The initial -stress algorithm [2,3] is used in the analysis. Three types of 
material stress-strain curves may be modeled: elastic-perfectly plastic, 
Ramberg-Osgood, and multilinear representations. Some unique features of 
the code are the computation of mixed-mode strain-energy release rates (G) 
for elastic solids using a 3D virtual -crack-closure technique (VCCT) [4], 
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the calculation of the J-integral for elastic-plastic materials using the 
equivalent domain integral (EDI) method [5-7], the capability to extend the 
crack under monotonic or cyclic loading, and the capability to close or open 
the crack faces during cyclic loading [8,9]. 

ELEMENT FORMULATION 

This section presents the element definition, which will be useful in 
data preparation, and the reduced-shear integration scheme [1] that is used 
to improve the bending performance of eight-noded elements. 

Isoparametric element formulation, given in many standard books (for 
example [2]), is used to define the shape functions and generate the element 
stiffness matrix. Figure 1 shows a 2-unit cube mapped into a general 

hexahedron element in the Cartesian coordinate system. The element shape 
function is determined for the cube and then transformed to any general 
element shape through coordinate transformation. Consider a local 
coordinate system £, r,, f for the cube as shown in Figure 1. The 

displacement or the coordinate at any point within the cube is defined by 
eight unknown q's corresponding to the eight nodes as 

<f> = aj + c*2 £ + + a^T} + a + <*g (1) 

Applying Equation 1 at all nodal points, a relation between the nodal 
coordinate values and is written in a matrix equation as 

<£ e = C a 

where the transpose of and a matrices are defined as 


( 2 ) 




{a}^ = (aj, a^, • • , Og) 

and C is a 8 by 8 square matrix consisting of coordinates (-1 or 1) of the 
eight nodes in the cube. Equation 2 is rewritten as 

a = C" 1 * e (3) 

Substituting Equation 3 into Equation 1 gives 

<£ = L a = L C 1 4> e (4) 

where 

t = [l, 17, £*7 > » 7 f » 

Thus, the element shape functions are defined by 

<t> = N f = [N p N 2 , . . , Ng] / (5) 

where 

N = L C" 1 (6) 

The element stiffness and consistent load vector evaluation follow the 
same procedures as given in Reference 1. In the element local coordinate 
system, the first two node numbers (1,2) in the element connectivity define 
the »?-axis and the second and third nodes (2,3) define the £-axis in the 
parent cell (cube). The coordinate system follows the right-hand rule. 
Therefore, the element face definition for surface loading (uniform pressure 
or stress) and corresponding program variable (IEFAC) is 


Face No. 

Planes 

Program 

Variable 

1 

i = -1 

IEFAC(l) 

2 

c = +1 

IEFACj 

2 ) 

3 

n = -1 

IEFAC I 

[3) 

4 

r; = +1 

IEFAC 1 

[4) 

5 

r = -1 

IEFAC ( 

5) 

6 

f = +1 

IEFAC(6) 
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As reported in References 1 and 2, the bending deformations of a linear 
element can be improved by performing reduced integration on shear strains. 
The procedure is straight forward for 2-D problems. However, for 3-D 
problems, a special procedure is needed to restrain rigid body deformation 
modes. The procedure outlined in Reference 1 is used in this program and is 
presented here for a two-point Gauss quadrature integration scheme. Each of 
the three shear strains («^, and e ^ ) are evaluated at two integration 

points as given in the following table, instead of at all eight Gauss 
points. 


SHEAR STRAIN INTEGRATION POINTS 

1 

**> 



o 

II 

II 

O 

H 

II 

C* 

1 

O 

II 

u 

ii 

r - -1/73 

t - -1/73 

v = -1/73 

f - +1/73 

€ - +1/73 

n = +1/73 


The integration points in the above table refer to the local coordinate 
system in the parent cube (or brick type) element. Care must be taken to 
avoid using highly skewed elements, especially when the reduced integration 
scheme is used. 


ELASTIC AND ELASTIC-PLASTIC ANALYSIS 
The analysis in the ZIP3D code consists of two parts. First, an 
elastic analysis is performed. Stresses, strains, nodal displacments, nodal 
forces and fracture parameters (G,J) are printed out as requested by the 
user. The program terminates if only an elastic analysis is requested. 
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If an elastic-plastic analysis is requested, the stresses, strains, nodal 
displacements and nodal forces for the previous elastic analysis are scaled 
to the incipient yield condition (that is, the highest stressed point is 
scaled to match the specified elastic limit of the material) and an 
incremental plasticity analysis is performed (see section on Solution 
Method). The selection of the elastic limit is not crucial in the analysis 
because the stress and strain state will always follow the input stress- 
strain curve (within a user specified tolerance). The elastic limit is used 
to establish the incipient yield condition and to obtain the incremental 
load (or incremental displacement) value. The finite-element formulations 
for the elastic and elastic-plastic analyses used in this program are given 
in Reference 1 and 2. An isoparametric finite-element formulation [2] was 
used to calculate the element stiffness matrix and consistent load vector 
due to distributed loads acting on the element faces. Element stiffness 
matrixes are assembled in a variable-band width stiffness matrix. A 
variable-band width linear equation solver given in Reference 8, based on 
Choleski's decomposition technique, is used to solve for the unknown nodal 
displacements. The solution algorithm and element stiffness generations are 
formulated differently from standard procedures to increase the vector 
lengths so that computations will be faster on vector processing computers. 

Figure 2 shows a flow chart for the elastic analysis in ZIP3D. The 
program has the option to impose symmetry boundary conditions on six 
different planes, as requested by the user. On the x = 0 and x = WIDTH 
planes, the displacement u = 0; on y = 0 and y = HEIGHT planes, the 
displacement v = 0; and on z = 0 and z = THICK planes, w = 0 can be 
specified. Also, nodes on the uncracked region of the crack plane are 
automatically restrained for FE models where the crack plane is on y = 0 
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Elastic-plastic analysis 


Figure 2.- Flow chart for elastic analysis 
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plane and the crack front is straight and normal to the x-axis, otherwise, 
the user must restrain the appropriate nodes to model the uncracked region. 

The parameter which controls whether an elastic or elastic-plastic 
analysis is performed is NEP (specifying NEP = 0 performs only an elastic 
analysis or NEP - n a positive number (n) performs an elastic-plastic 
analysis. During a plastic analysis, stresses, strains, displacements and 
nodal forces, as specified by the user, are printed out after every NEP 

(n^) load or displacement steps. 

The elastic-plastic analysis is based on incremental plasticity theory 
using the von Mises yield criterion, isotropic hardening, and Drucker's flow 
rule. The initial -stress algorithm developed by Zienkiewicz et.al. [3] is 
used in the analysis to satisfy the equations of equilibrium and the stress- 
strain relation. The von Mises yield criterion for a 3-D stress state is 
defined by: 

FU) = - a (7) 

where 

a eff “ «('xx " a yy^ + ^ a yy ’ a zz^ + ^ a zz ” a xx^ 

t 6, 2 xy + 6j yz ♦ 6„2 x ]/2) 1 / 2 . 

The stresses a . ■ represent the 3-D stress components at any Gauss point and 

' J 

a is the current flow stress at the same Gauss point. If F(<7) < 0, the 
material is in an elastic stress state and if F(a) > 0, the material is in 
a plastic state. A radial -return technique [10] is used to bring the stress 
field to yield surface, F(cr) = 0. The difference between the total 
incremental stress and the true incremental stress gives the residual 
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stresses which contribute to the plastic load vector Q. A flow chart for 
the step-by-step procedures in the elastic-plastic analysis is presented in 
the Figure 3. Convergence of the solution is defined by 

(a eff - a)/a < ERIT (8) 

where ERIT is the user specified accuracy. The elastic-plastic analysis is 
continued by applying incremental load (or displacement) steps, as specified 
in the input data by PCT, until the desired load (or displacement) is 
reached. PCT is a percentage of the load (or displacement) that causes 
incipient yielding on the highest stressed element. At the end of a 
specified number of load steps (NEP), stresses, strains, displacements, 
forces, and J-integrals are calculated and printed out. 

Three types of material stress-strain curves can be used in this 
analysis: elastic-perfectly plastic, Ramberg -Osgood, and multilinear 
representations. These three types of stress-strain representations are 
shown in Figure 4. 

The ZIP3D program can be used to analyze stationary cracks as well as 
growing cracks. In stationary crack problems, the analysis can be continued 
for any type of loading and unloading cycles by inserting a series of load 
factor (ratio of maximum load to initial applied load) input lines. The 
program stops when it encounters "HALT" in columns 12-15 in the load-factor 
input line. 

SOLUTION METHOD 

The application of the finite-element method to problems involving 
linearly elastic materials is straightforward because the material 
properties are constant and only one solution is required to obtain 
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Ramberg-Osgood 
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displacements for the elastic structure. However, for elastic-plastic 
structures the coefficients in the stiffness matrix are functions of 
loading. Thus, the displacements usually are obtained by applying small 
load increments to the structure and either updating the coefficients of the 
stiffness matrix or applying an "effective" plastic-load vector after each 
load increment (initial -stress method). The latter technique was used in 
the ZIP3D code. 

In general, the matrix equation which governs the response of a 
discretized structure under loads which cause plastic deformations is 


[K e ] {U}j = (P) 1 + (Q}j:j (9) 

where [K g ] is the elastic stiffness matrix, (U) is the generalized nodal 

displacement vector, (P) is the applied load vector, and {Q} is the plastic 
load vector. The superscript i in Equation 9 denotes the current load 
increment, and i-1 denotes the preceding increment. After each load 
increment an iterative process is required to stabilize the plastic-load 
vector. The subscript I in Equation 9 denotes the current iteration, and 
1-1 denotes the preceding iteration. The elastic element stiffness 
matrixes are assembled in a variable-band width storage routine. A 
variable-band width linear equation solver given in Reference 8, based on 
Choleski's decomposition technique, is used for the solution. The solution 
algorithm and element stiffness generations are formulated differently from 
standard procedures to increase the vector lengths so that the computations 
will be faster on vector processing computers. In the initial -stress 
method, the solution to an elastic-plastic problem is obtained by applying a 
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series of small load increments to the structure until the desired load is 

reached, {P}^ = {P} i_1 + {dP}. The incremental load vector {dP} is chosen 
as a percentage of the applied load (or applied displacement) that causes 

the highest stressed element to yield. During the i** 1 load increment a 
purely elastic problem is solved, and the increments in total strain {de} 
and corresponding elastic stress {da g } are computed from the displacements 

for every element. Because of material nonlinearity the stress increments 
are not, in general, correct. If the correct stress increment for the 
corresponding strain increment is (da), then a set of body forces or 
plastic-load increment {dQ} caused by the ''initial" stress {da 0 } (= {da g } - 

{da}) is required to maintain the stress components on the yield surface. 
The plastic-load increments on each element are computed from 


{dQ} = J [B] T {da“} dV 


( 10 ) 


where [B] is the strain-displacement matrix and V is the volume of the 
element. For elements which are in an elastic state or unloading from a 
plastic state, {dQ} = 0. The total plastic-load vector is computed as 


<Q}j = {Q)j;l + (dQ) 


( 11 ) 


The new force system {Q}J is added to the applied load vector in Equation 

(9) and a new set of displacements is obtained. The iteration process is 
repeated until the stress state is sufficiently close to the stress-strain 
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curve of the material (see Eqn. 8). Usually, 5 to 15 iterations are 
required to stabilize the plastic-load vector. However, for configurations 
which have large strain gradients, more iterations are required. In order 
to reduce the number of iterations, a relaxation technique [2] was 
incorporated into the plastic analysis program by using 

(Q}j - (Q)j:J + 9 {dQ} (12) 

where g is the relaxation parameter (RP). Because the displacements from 
the preceding increment or iteration are used to compute the plastic-load 
increment, the plastic-load vector is underestimated. Thus, the relaxation 
parameter is used to increase the plastic-load vector and, consequently, 
increase the rate of convergence. A value of RP = 1.5 has been found to 
give accurate solutions in about two-thirds of the CPU time as that for RP = 
1. Higher values of RP are not recommended. 

CRACK EXTENSION METHODS 

Two types of crack extension criteria are used to simulate crack front 
movement. They are the critical crack-tip-opening-displacement criterion 
and a specified load criterion. In the first criterion, the crack is 
extended when the crack-tip opening (v) displacement (CTOD) at the node 
immediately behind the crack-tip node exceeds a specified value (SCRIT). In 
the critical CTOD approach, the CTOD value will depend upon element size. 
Therefore, the element size and pattern around the crack front must be the 
same for other crack sizes or other crack configurations. In the second 
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criterion, when the applied load equals or exceeds a specified input value 
the crack is extended. 

The program has the option to extend the crack either to create a new 
crack surface during loading or to close and open the crack surfaces during 
cyclic loading. Crack extension is restricted to symmetric problems, where 
the crack plane lies on the y = 0 plane and the crack front is straight. 
Symmetry boundary conditions on the y * 0 plane are enforced by using 
"rigid” springs on the uncracked region (springs with a very large 
stiffness). The spring stiffness, K $ , is added to the diagonal stiffness 

coefficient in [K g ]. During crack extension, the crack front extends by one 

element size and the crack front remains straight after extension. FE 
models multilayered in the z-direction, such as that shown in Figure 5(a), 
are recommended for crack extension analyses. The crack extension process 
is shown in Figure 5(b). The program identifies nodes on the crack plane 
and those in the immediate region around the crack front. To extend the 
crack, the transverse (y-direction) restraint at each of the nodes on the 
crack front is replaced by a nodal force (F) acting in the y-direction; the 
factorized stiffness matrix is modified to account for the change in the 
boundary condition associated with each of the separated crack front nodes. 
This keeps the new crack faces closed and the stresses and displacements in 
the body are identical to the previous solution. The equations of 
equilibrium are then given by 

[K43 (U)J - (P) 1 + (Q}j:J + (F) (13) 
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CRACK 


Crack front 


(a) Finite-element idealization at the crack plane 



| — ► Crack extension 
After 


(b) Crack extension process 




(c) Crack closure process 


Figure 5.- Crack extension and closure process for a straight crack front problem 
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where [K^] = [K g ] + [K $ ]. [K $ ] Is the diagonal matrix for intact springs. 

The change in boundary conditions associated with crack advance is achieved 
by replacing the rigid spring with a zero-stiffness spring. See Reference 8 
for details on modifing the Cholesky factors in the stiffness matrix. The 
crack-front nodal forces {F} are now released incrementally in steps 
(usually 4 to 5 steps, a value NLM specified by the user in the input data). 
After each incremental nodal -force (dF) change, a plasticity analysis is 
performed. When the nodal forces on the crack-front nodes become zero ({F} 
» 0 in Eqn. 13), a new crack surface equivalent to one-element length is 
created and the crack front has advanced by the same amount. In this 
process, residul plasticity deformations are left behind as the crack 
advances. 

To model crack closure during unloading, the crack-surface nodal 
displacements are monitored and if any nodes come in contact or overlap, a 
rigid spring is inserted to carry the compressive loads (see Fig. 5(c)). 
The equations of equilibrium are 


[KJ] (U)j - (P) 1 + {Q}j;J (14) 
where [K^] - [K'] + [y. 

The crack opening and closure operations are performed in subroutines 
BREAK and CONTACT, respectively. Although the program is specialized for 
straight crack fronts, the program can be modified to analyze extension of 
curved crack fronts. However, the growth criterion would have to be 
modified and would now be a function of location through the plate 
thickness. 
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EVALUATION OF FRACTURE PARAMETERS 

The most commonly used fracture parameters are the stress- intensity 
factor and strain-energy-release rate for linear-elastic fracture mechanics 
and the J-integral for nonlinear fracture mechanics. Two methods of 
evaluating these parameters for cracks in solids are presented: the virtual - 
crack-closure technique (VCCT) [4] and the equivalent-domain-integral (EDI) 
method [5-7]. The VCCT is used for elastic mixed-mode fracture problems and 
the EDI method is used for both elastic and elastic-plastic fracture 
problems. The EDI method used here calculates only the total J-integral and 
not its components for mixed-mode crack problems. The two methods are 
described in the following sections. 

Strain-Energy-Release Rates and Stress-Intensity Factors 

The three-dimensional virtual -crack-closure technique [4] used here is 
an extension of the 2-D technique of Rybicki and Kanninen [11]. The present 
method accounts for the variation of the stress-strain field along the crack 
front. The method calculates the energy required to close the crack 
surfaces on a segment of the crack front. Figure 6 shows the cross-section 
and plan view of a FE idealization near the crack front. The crack front is 
divided into several segments as shown in Figure 6. A segment is the 

thickness of an element. Consider the i** 1 segment along the crack front 
defined by nodes p and q. The Cartesian coordinate system (x, y, and z) and 
the corresponding displacements (u, v, and w) are defined as shown. The 
local crack front coordinate system is represented by Xj, and x^. The 

FE idealization around the crack front should be nearly orthogonal and the 
element pattern and sizes immediately behind and ahead of the crack front 


y,v 



(b) Top crack face 


Figure 6.- Finite-element idealization at a segment of the crack front 
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should be nearly identical. The elements, for example, e^, e 2 > e 3 , and e^, 

surround the crack front. Nodes p and q are on the crack front. Nodes m 
and n are behind the crack front on the top face; and nodes m' and n' are on 
the bottom crack face. At each of the crack-front nodes there is a self- 
equil ibriating force system as shown in the insert in Figure 6. If the 
element size around the crack front is small compared to the crack length 
(ratio less than or equal to 1/10), the total strain-energy-release rate 
(Gj) is given by 

G T = [F V n ( U m _ u m /) + F wri (v m - v m ,) 

I L xp v m m 7 yp v m m ' 

* F zp ■ V> + F xq < u n - V> < 15 > 

+ F yq < v n ‘ V> + F zq < w n ’ 

where u, v, and w are displacements in x, y, and z directions, respectively. 

F • , F ., and F • refer to nodal forces acting in the x, y, and z directions 
aj y j l j 

respectively. The subscript j refers to nodes p and q. Lengths Aj and a 3 

are the length and width of the crack front element in the x^ and x 3 

directions, respectively. Gy calculated by the Equation 15 is an average 

value over the crack front segment (p to q) and, hence, it is assumed to be 
valid at the mid-point of the segment. 

Components of strain-energy-release rates in the opening (Gj), 

shearing (Gjj)> and tearing (Gjjj) modes can be evaluated by resolving the 

nodal forces and displacements into components in the local coordinate 
system (x^, X£, and x 3 ). Applying the virtual crack closure principles to 
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components of forces and displacements for opening, shearing, and tearing 


modes, the Gj, Gjj, and Gjjj equations are written as 


G I * * F 2p ( v 2m ' v 2m') + F 2q ^ v 2n " v 2n' )>/( 2 a i a 3 } (16) 

G II * ^ F lp ( u 2m * u 2m^ + F lq ( u 2n ' u 2n' )V(2 A i A 3 ) (17) 

G III = fF 3p ^ w 2m ' w 2m') + F 3q ( w 2n ' W 2n'^/^ 2A 1 A 3^ ^ 18) 

Corresponding stress-intensity factors in modes I, II, and III are given by 

Kj = 7“Gp* (19) 

K II = (20) 

K III * J G III E/(1 + (21) 

where 

2 

E* = E/(l-i/ ) for plane-strain condition 
E* » E for plane-stress condition 
E = Young's modulus 
v = Poisson's ratio 


J-Integrals 

The Equivalent Domain Integral (EDI) method is used in ZIP3D to obtain 
the total J-integral for cracked elastic and elastic-plastic bodies. The 
total J-integral at any point along the crack front in a 3D cracked body is 
defined as an integral over a closed surface around the crack front [12] 
(the tube represented by the broken line in Fig. 7) as 
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J = 


Lim 1 
r/A -* 0 - 

A -» 0 A 


[w nj - au^/axj nj] dA 


( 22 ) 


A r + o 1+ o 2 


where w is the stress-work density, defined as 




dc ij 


(23) 


In Equation 22, a. . and e.. are stress and strain tensors on the surface of 

I J • J 

the tube, u i is the displacement vector, n- is the j th component of the unit 

normal vector on the surface, a is the projected length of the crack front 
along the axis, and r is the radius of the tube over which the integral 

is evaluated. The indices i and j take the values 1, 2, and 3. Recently, 
the surface integral in Equation 22 was modified to a volume integral, 
called the equivalent domain integral [5-7], for ease of implementation in a 
finite-element analysis and accurate evaluation of the integral. This was 
done using Green's divergence theorem and de Lorenzi's s-function [13]. For 
traction-free crack faces, Equation 22 is written as a volume integral 
between the inner tube A r (broken line) and any other closed tube (solid 

line) enclosing A r (see Fig. 7) as 


j . - i 

f 


[w as/axj - auj/axj as/axj] dv 
[ aw/axj - tr.j Se^/axj ] s dv j 


(24) 
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Figure 7.- Domain around the crack front 


where s is any arbitrary but continuous function with a characteristic value 
of one on the inner tube and zero on all outer surfaces, and f is the 
integrated value of s along the x^- axis. Equation 24 can be used for 

isotropic or anisotropic as well as linear or nonlinear materials. In the 
present incremental elastic-plastic analysis, the stress-work density, 
stresses, and strains were evaluated at each load step. Then, the J- 
integral was evaluated at specified load steps from the accumulated stress- 
work density, stress, total strain, and displacement fields. 

In the ZIP3D code, the user inputs only the element numbers within a 
selected domain (between A and A r surfaces and within the volume V) and the 

node numbers on the "inner" surface A r at which s = 1. The program 

calculates s- values for all other nodes and automatically sets up the 
appropriate s-functions and f-integrals (Eqn. 24) for each domain specified. 
Because the 8-node hexahedral element is a linear displacement element, only 
linear s-functions are used. Details of general s-functions are given in 
References 5 and 7. 

The selection of the domains for evaluating the J-integrals depends 
upon several factors for 3-D crack configurations. The important factors 
are the gradient of the strength of the singularity along the crack front, 
radii of the inner and outer surfaces of the domain, and the finite-element 
modeling in the crack front region. Some simple guidelines are stated here. 
Further details are given in Reference 7. The crack front need not to be 
modelled with singularity elements and the FE idealization at the crack 
front can be either rectilinear or polar. The radius of the inner surface 
of the domain should be small (less than 1/10 of the crack length) or zero. 
The radius of the outer surface of the domain should be less than 1/5 of the 
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crack length. Because the amount of input data increases with the number of 
elements in the domain, it is convenient to use one or two rings of elements 
around the crack front. If the FE mesh around the crack front is a polar 
mesh, the first ring of elements should give an accurate solution. However, 
for rectilinear meshes the second ring of elements is recommended. For an 
elastic-plastic analysis the second ring of elements is recommended. 
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NOMENCLATURE 


a surface crack depth 

B plate thickness 

[B] strain-displacement relation 

c crack length 

{dP} incremental applied load vector 

{dQ} incremental plastic-load vector 

(de) incremental strains 

(da e ) incremental elastic stress increments 

{da 0 } incremental initial stress 

E Young's modulus 

{F} crack front nodal force vector 

F(ct) von Mises yield criterion 

f area under s-function curve 

G, Gj strain-energy release rate 

G. mode i strain-energy release rate (i = I, II, III) 

g relaxation parameter 

H plate height 

J J-integral 

K. mode i stress-intensity factor (i = I, II, III) 

[K e ] elastic stiffness matrix 

[Kg] modified elastic stiffness matrix 

[K $ ] spring stiffness diagonal matrix 

m Ramberg -Osgood stress-strain curve power 

(P) applied load vector 
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{Q} 

s 


plastic-load vector 
applied stress 
s deLorenzi's s-function 

W plate width 

w stress-work density 

{11} generalized displacements 

u,v,w displacements in x-,y- and z-directions, respectively 
V element volume 

x,y,z Cartesian (global) coordinate system 

Xj,x 2 ,x 3 local crack front coordinate system 

a ^ constants used in defining element shape function 

a element length along crack front 

<f> element shape function 

c • ^ strain tensor 

k Ramberg-Osgood stress-strain coefficient 

v Poisson's ratio 

stress tensor 

a flow stress 

a eff effective stress 

a Q elastic (proportional) limit stress 

^ , » 7 , 5 * element local coordinate system 
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III. USER MANUAL FOR ZIP3D 

INTRODUCTION 

A line-by-line explanation of the input data for the ZIP3D code is 
presented herein. Each of the parameters used in the input data file is 
defined and explained. Input data preparation is explained with reference 
to a single-edge-crack-tension specimen, for the purpose of clarity. 
Because the finite-element (FE) idealization used in this model is coarse, 
the user is urged to use this data only as a guideline and not for producing 
an accurate solution. 

Figure 8(a) shows a single-edge-crack-tension specimen with crack 
length a, width W, thickness B, and height H. For symmetric loading, only 
one-fourth of the specimen is modeled in the FE analysis. Figure 8(b) shows 
the FE idealization at z = 0 plane and Figure 8(c) shows the 3-D 
idealization. The overall size of the FE model is defined by the WIDTH, 
THICK, HEIGHT parameter, see Fig. 8(c). The finite-element idealization at 
the crack plane and around the crack front are shown in Figures 8(d) and 
8(e), respectively. 

For ease of handling, the input data is specified in two files: an 
analysis file and a mesh file. The analysis file consists of the geometric 
definition of the cracked body, number of nodes and elements in the finite- 
element idealization, material data, boundary conditions, loading, G and J 
calculation data, plastic analysis data, and output specifications. Free 
format read statements are used in the program, except as noted. 

The mesh file consists of nodal coordinates (x, y, and z), element 
connectivity, and information on pressure loadings on each of the element 
faces. An option for reduced-shear integration or full integration for each 
of the elements is specified by the variable INDXEL. Because it is 


30 


|y 



(a) Single-edge-crack 
tension specimen 


tr 


2 3 4 





r— T| - ~ 

U 

□ 

UJ 

E 

4 


6 

7 

8 

9 


E 

11 

6 

12 

B 

13 

8 

14 

k 


im 

m 

m 

MB 

■ 

/ 

Crack tip — ' 
(b) Mesh at z 

V * w 

= 0 plane 


10 

15 

20 


t y 



Face 

ISYMPL 

o 

II 

X 

1 

y = o 

2 

z = 0 

3 

x = WIDTH 

4 

y = HEIGHT 

5 

z = THICK 

6 


(c) Three-dimensional idealization with element and node numbers 



Crack front 



x 


(d) Crack plane idealization (e) F-E idealization at the mid-plane 


Figure 8. Example problem and finite-element idealization. 
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useful to keep the original node numbers as they are input and also to print 
out the results in the original numbering scheme, an option to input an 
optimized nodal numbering scheme is provided. The user must provide the 
optimized nodal numbering scheme. Several optimizing codes are available in 
the literature and are not discussed here. An optimized numbering scheme 
will reduce the storage and solution time. However, output data will be 
printed out in the original nodal numbering scheme. A line-by-line 

description of analysis and mesh files are presented in the following 
sections. 

PROGRAM SIZE CONTROL 

The code was written in ANSI standard Fortran 77 language for vector 
processing computers. The code was originally developed for the Cyber 205 
super-computer and has been converted to CRAY2 computers. The program 
utilizes several CRAY2 intrinsic functions to speed up the computation. The 
program can be implemented on any other super-computer or workstation by 
appropriately changing the intrinsic functions. The complete program is 
written using variable dimensions and parameter statements as 
PARAMETER (LNSTIF = 11000000) 

PARAMETER MXNOD = 6000, MXNEL = 5000, MXBND = 3*MXNOD) 

PARAMETER (MXPN0D = 2000, MXBCND = 975, MXNL1 = 10, MXMAT = 3) 

The program size is controlled by the eight variables defined in these 

parameter statements. They are: 

LNSTIF - maximum profile length of stiffness matrix (11,000,000) 

MXNOD - maximum number of nodes (6000) 

MXNEL - maximum number of elements (5000) 

MXBND - maximum bandwidth of assembled stiffness matrix (3*MXN0D) 

MXPNOD - maximum number of loaded nodes (2000) 

MXBCND - maximum number of restrained nodes or maximum number of 
specified displacement degrees-of-freedom (975) 
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MXNLI - maximum number of layers in thickness (z-direction) + 1 (10) 
(only used in crack-extension analysis) 

MXMAT - maximum number of materials (3) 

The numbers within the parenthesis represent the initial values set in 
the program. Whenever the input data exceeds any of these specified values, 
the program flags the exceeded variable and stops. The program size can be 
increased by changing the value of any variable in all parameter statements 
in the program. 


ANALYSIS (INPUT DATA) FILE 


Line No. Parameter Format 

1. TITLE - Analysis Title (80-characters) 20A4 

2. MFILE - Mesh Filename A10 

3. CRACK, WIDTH, THICK, HEIGHT, DAX, SCALE, PRSUR 

CRACK - crack length, applicable for problems 
with crack plane on the y = 0 plane and 
the crack front is normal to the global 
x-axis of the specimen, see Fig. 8(c). 

(set to zero for uncracked body) 

WIDTH - width of the FE model, measured parallel 
to the global x-direction 

THICK - thickness of the FE model measured parallel 
to the global z-direction 

HEIGHT - height of the FE model, measured parallel 
to the global y-direction 

DAX - element length in x-direction at crack front, 
applicable for crack extension analyses only 

SCALE - scale factor, globally scales model to 
desired size 

PRSUR - intensity of pressure applied on any element 
face specified in the mesh file 
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PRINT AND PROGRAM OPTIONS 


4. LPRIT, NNODE, NELM, N LAYER, NEP, IOPTON, KVCCT, JEDI 

LPRIT = 0 - short output 

= 1 - detailed output 

NNODE - number of nodes in the model 

NELM - number of elements in the model 

NLAYER - number of layers in z-direction, 

layers must be repetitions of first 
layer, see example in Figure 8. 

(used only for controlling output) 

NEP -number of load increments between output 
of results in the plastic analysis. 

NEP = 0 - elastic analysis of cracked body, with 
crack plane on y = 0 plane and the 
crack front is normal to the x-axis 

NEP * -1 - elastic analysis of general solid with 
or without a crack (crack plane not 
on y = 0 plane) 

NEP = n - elastic-plastic analysis of cracked body 

with output at every n^-load increment 
(crack plane on y = 0 plane) 

NEP = -n - elastic-plastic analysis of an uncracked body 
or crack plane is not on y = 0 plane 
with output at every (n-1) load increments 

IOPTON = 0 - stress output at element centroid 

IOPTON = 1 - stress output at Gauss points 

IOPTON = 2 - stress output at nodal points 

KVCCT = 1/0 - Yes/No calculate G using the 3D VCCT 
method 

JEDI = 1/0 - Yes/No calculate J using the 3D EDI 
method 


MATERIAL PROPERTY DATA 
5. 


6 . 


NMAT - number of material types 
YOUNG, POIS, SIGYS, AM, ROM 
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YOUNG - elastic modulus 
POIS - Poisson's ratio 

SIGYS - elastic (proportional) limit stress 

AM - specifies type of material stress-strain 
representation 

AM = 0 - elastic-perfectly plastic 

AM > 0 - Ramberg -Osgood exponent 

AM = -1 - piecewise linear representation 

ROM - Ramberg-Osgood constant as defined in 
the stress-strain equation 

f a T AM 

£ = YOUNG + l ROM ' 

If AM = -1 continue, otherwise go to 9. 

7. NSEGMT - number of piecewise linear segments in stress- 

strain curve description 

8. YSTRS(i), YSTRN(i) - stress and total strain at end of 

segment i on the stress-strain 
curve (see Fig. 4) 

[Repeat line 8 NSEGMT times] 

9. NSET, (NBEGIN(i), NEND(i) , NINC(i), i = 1 to NSET) 

NSET - number of sets of elements having material 
properties specified on lines 6-8. 

NBEGIN(i), NEND(i), and NINC(i) - Beginning, ending, and 

increment in element number 
for material property set i 


[Repeat lines 6 to 9 NMAT times] 

SYMMETRY BOUNDARY CONDITIONS 

10. NSYMPL - number of symmetric conditions 

(only used when all nodes on 
a plane are to be restrained; 
otherwise set NSYMPL = 0) 

11. (ISYMPL(i), i = 1, NSYMPL) (skip if NSYMPL = 0) 

- plane numbers (see Fig. 8(c)) 

= 1 x = 0 

= 2 y = 0 
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= 3 z = 0 

=4 x = WIDTH 

=5 y = HEIGHT 

=6 z = THICK 

FIXED. DISPLACED AND LOADED NODES 

12. NFIX, NLOAD, NSPD 

NFIX - number of restrained nodes 

NLOAD - number of nodes with specified loads 

NSPD - number of nodes with specified displacements 

13. JFIX, MU, MV, MW (If NFIX > 0 continue, otherwise go to 14) 
JFIX - restrained node number 

MU = 1/0 - restrained/unrestrained in x-direction 
MV = 1/0 - restrained/unrestrained in y-direction 
MW = 1/0 - restrained/unrestrained in z-direction 
[Repeat line 13 NFIX times] 

14. NODS, K, DISP (If NSPD > 0 continue, otherwise go to 15) 

NODS - node number of the specified displacement 

K = 1 - displacement in x-direction 

K = 2 - displacement in y-direction 

K = 3 - displacement in z-direction 

DISP - specified displacement 
[Repeat line 14 NSPD times] 

15. NODLOD, PX, PY, PZ (If NLOAD > 0 continue, otherwise go to 16) 
NODLOD - node number 

PX - load in x-direction 

PY - load in y-direction 

PZ - load in z-direction 

[Repeat line 15 NLOAD times] 
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CRACK EXTENSION OPTIONS 


16. NTYP, NLM, SCRIT, RP, ACURCY 

NTYP - specifies crack extension criteria 

NTYP - 0 - crack extension at specified load (see line 30) 

NTYP * 1 - crack extension at specified CTOD 

NLM - number of steps to release crack-tip force (see Eqn. 13) 

SCRIT - critical value of CTOD for crack extension 

RP - relaxation parameter to accelerate convergence 
of elastic-plastic analysis: specify value of 

RP between 1 and 2 (RP = 1 is normal analysis 
without acceleration.) 

ACURCY tole^nce; between SCRIT and calculated crack- 
tip opening displacement to extend crack 

NOTE: NTYP, NLM, SCRIT, and ACURCY are required only for crack 

extension analysis (use dummy values for other analyses) 


PI ASTICITY OPTIONS AND NODAL/ELEMENT OUTPUT INFORMATION 
17. PCT, ERIT, MAXIT, N0DE1, NODE 2 , NELE1, NELE2 

PCT - increment load as percentage of yield load 
(aP = PCT * load at elastic limit) 

ERIT - allowable error for stress convergence 
criteria (terminates plastic-load vector 
iterations), see Eqn. 8. 

MAXIT - maximum number of iterations allowed 
(stops divergent solutions) 

NODE1 and N0DE2 - beginning and ending node numbers for 

output of nodal displacements and forces. 

Nodal stresses are also output if I0PT0N = 2. 

Set N0DE1 = NODE2 = 0 to print out only fracture mechanics 
parameters (nodal displacements, forces and stresses 
will not be printed) 

NELE1 and NELE2 - beginning and ending element numbers for stress 

output (IOPTON = 0 or 1). If N LAYER > 1, 

NELE1 and NELE2 correspond to element numbers 
of the first layer. Stress output for 
corresponding elements in other layers 
is automatic. 
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G-CALCULATION BY VCCT 


If KVCCT = 0 Skip lines 18 through 21 

18. NGCAL - number of crack front segments were G is desired 

(a segment is the width of an element) 

19. JSYM - symmetry on/off switch 

JSYM - 0 - general unsymmetric crack problem 

JSYM = 1 - crack plane coincides with x = 0 plane of the model 

JSYM = 2 - crack plane coincides with y = 0 plane of the model 

JSYM = 3 - crack plane coincides with z = 0 plane of the model 

20. NGEL, ( LEG ( i ) , i = 1 to NGEL) 

NGEL - number of elements around crack front; top half of 
model for symmetric problems and all elements 
around the crack front for unsymmetric problems 

LEG(i) - element number 

21. (NGV(j) , j = 1 to NODV) , NGF(l) , NGF (2) 

NGV(j) - node numbers for nodes behind crack front segment. 

For symmetric problems, use the two nodes on the 
top crack surface (NODV - 2) and for unsymmetric 
problems, use both top and bottom crack surface 
nodes (NODV = 4) in the same order as specified 
in Figure 6 (m-n and m'-n'). 

NGF(k) - node number on the crack front segment (see Fig. 6, 
nodes p-q; k=l for node p and k=2 for node q). 

[Repeat lines 20 to 21 NGCAL times] 


J-CALCULATION BY EDI METHOD 

If JEDI = 0 Skip lines 22 through 29 

22. NEDIS - number of domains for EDI calculation of J-integral 

23. IROTYP - definition of local crack front coordinate system 

IROTYP = 0 - Cartesian coordinate system defined by 
three points 

IROTYP = 1 - polar coordinate system, where crack plane is 
on x-y plane and the angle (THETA) is measured 
from the x-axis. 
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24. x Q , y Q , z Q , Xj, y p z p x 2 , y 2 , z 2 , if IROTYP = 0 

coordinates of the three points 0, 1, and 2 (see Fig. 9(a)). 
Point 0 is on the crack front where J is evaluated. Point 1 
is on the crack plane and the vector from points 0 to 1 is 
normal to the crack front. Point 2 is on a plane normal to 
the crack plane and the vector from points 0 to 2 is normal 
to the vector from points 0 to 1. 

24. THETA, if IROTYP = 1 

angle between the x-axis and a vector drawn normal to the 
crack front, measured on the crack plane, from a point at 
which the J-integral is required (see Fig. 9(b)). 

25. NELEDI, LAYER 

NELEDI - number of elements in the domain 

LAYER - number of crack front segments in the domain 

26. (ELNEDI(i), i = 1 to NELEDI) 

ELNEDI(i) - element numbers for all elements within domain 

27. NODEDI - number of nodes in the domain with s = 1 (s-function) 

28. (SNODS(i), i = 1 to NODEDI) 

SNODS - node numbers with s = 1 in the domain 

29. [SNODX(i), i = 1 to (LAYER + 1)] 

SNODX - node numbers on the crack front element 
within a domain 


APPLIED LOAD/DI SPLACEMFNT FACTOR AND TERMINATION OPTION 

30. P, WORD F10.4, 1X,H4 

P - desired load factor 

WORD = GROW - extend crack at P 

= "blank space" - contines analysis with 
no crack growth 

[Repeat line 30 to describe monotonic or cyclic load-factor history; 
such as Pj, P 2 , P 3 and etc.] 

31. WORD 1 IX , H4 

WORD = HALT - terminates the analysis 


39 



40 


FINITE-ELEMENT MESH FILE 


Line No. Parameter Format 

1 . JRENUM 

JRENUM = 0 - original node numbering scheme is 
used in the analysis 

JRENUM = 1 - renumbering string for mesh optimization 
is provided 

2. (JNEW(i), i = 1 to NNODE) if JRENUM = 1 

JNEW(i) - string of optimized node numbers 

3. NODE, x, y, z 

NODE - node number 

x, y, z - Cartesian coordinates of the node NODE 
[Repeat line 3 NNODE times] 

4. IE, [MODE ( I E , j ) , j=l to 8], INDXEL, [IFACE(j) , j=l to 6] 1615 

IE - element number 

MODE - nodal connectivity defined in a counter-clockwise 
direction using the right-hand rule (1, 2, 3, 4, 

5, 6, 7, 8 as shown in Fig. 1) 

INDXEL = 0 - reduced shear integration, select this option 
for bending dominant problems 

INDXEL = 1 - full integration, select this option 
for tension dominant problems 

IFACE(j) = 0 - jth face is not pressure loaded 

I FACE ( j ) = 1 - jth face is uniformly loaded with a pressure of 
intensity PRSUR (Line #3 in Analysis File). 
Element face numbers are defined in Figure 1 
with reference to the local coordinate system. 

[Repeat line 4 NELM times] 
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EXAMPLE DATA FILES 


The following analysis and mesh files are for an elastic analysis of a 
single-edge-crack-tension specimen subjected to uniform stress at the ends 
of the specimen. Only one-fourth of the specimen is idealized by 8-node 
elements. The specimen configuration and details of the finite-element mesh 
are shown in Figure 8. This example is for a coarse mesh, hence it should 
be used as a guideline only, not for generating accurate solutions. 


ANALYSIS FILE 


EXAMPLE DATA FOR ELASTIC ANALYSIS OF SINGLE-EDGE-CRACK-TENSION SPECIMEN 
MESHFILE 

.5, 1., .5, 1., .2, 1. 

0, 60, 24, 2, 0, 1, 1 
1 

1., .3, .5, -1., 0.0 
3 


0.5, 0.5 
0.7, 0.75 

0. 7, 1.00 

1, 1, 24, 1 
1 

3 

1 , 0 , 0 , 1.0 
1 , 1 , 0 , 0 
0, 1, 1., 1., 0.99 
1., .01, 100, 1, 60, 1, 12 
2 $ G by 

0 


$ number of symmetric planes 
$ z = 0 plane 

$ node 1 u=Q 

$ not applicable for elastic analysis 

$ two layer analysis 

VCCT method at 2 crack front segments 


2 , 10 , 11 

17, 37, 18, 38 
2, 22, 23 

37, 57, 38, 58 

3 

0 

0 ., 0 ., 0 ., 1 ., 0 ., 
2 , 1 
10, 11 
1 

18 

18, 38 
0 

0 ., 0 ., 0 ., 1 ., 0 ., 
4 , 2 

10, 11, 22, 23 
1 
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$ set # 1 
$ set # 2 

$ J by EDI method at 3 locations on the crack front 
$ set # 1 at z=0 
0 ., 0 ., 1 ., 0 . 


$ set #2 at z= THICK/2 

0 ., 0 ., 1 ., 0 . 
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18, 38, 58 

0 $ set #3 at z=THICK 

0 ., 0 ., 0 ., 1 . , 0 . , 0 . , 0 . , 1 . , 0 . 

2 , 1 

22, 23 

1 

58 

38, 58 

1.0 HALT $ required only for elastic-plastic analysis 


MESHFILE 


1 

0 . 

1 . 

0 . 

2 

.3 

1 . 

0 . 

3 

.5 

1 . 

0 . 

4 

.7 

1 . 

0 . 

5 

1 . 

1 . 

0 . 


$ x, y and z coordinates of nodes 


— (continue) 

56 0. 0. 1. 

57 .3 0. 1. 

58 .5 0 . 1 . 

59 .7 0. 1. 

60 1. 0. 1. 


1 

1 

6 

7 

2 

21 

26 

27 

22 

1 

2 

2 

7 

8 

3 

22 

27 

28 

23 

1 

3 

3 

8 

9 

4 

23 

28 

29 

24 

1 

4 

4 

9 

10 

5 

24 

29 

30 

25 

1 

5 

6 

11 

12 

7 

26 

31 

32 

27 

1 

6 

7 

12 

13 

8 

27 

32 

33 

28 

1 

7 

8 

13 

14 

9 

28 

33 

34 

29 

1 

8 

9 

14 

15 

10 

29 

34 

35 

30 

1 

9 

11 

16 

17 

12 

31 

36 

37 

32 

1 

10 

12 

17 

18 

13 

32 

37 

38 

33 

1 

11 

13 

18 

19 

14 

33 

38 

39 

34 

1 

12 

14 

19 

20 

15 

34 

39 

40 

35 

1 

13 

21 

26 

27 

22 

41 

46 

47 

42 

1 

14 

22 

27 

28 

23 

42 

47 

48 

43 

1 

15 

23 

28 

29 

24 

43 

48 

49 

44 

1 

16 

24 

29 

30 

25 

44 

49 

50 

45 

1 

17 

26 

31 

32 

27 

46 

51 

52 

47 

1 

18 

27 

32 

33 

28 

47 

52 

53 

48 

1 

19 

28 

33 

34 

29 

48 

53 

54 

49 

1 

20 

29 

34 

35 

30 

49 

54 

55 

50 

1 

21 

31 

36 

37 

32 

51 

56 

57 

52 

1 

22 

32 

37 

38 

33 

52 

57 

58 

53 

1 

23 

33 

38 

39 

34 

53 

58 

59 

54 

1 

24 

34 

39 

40 

35 

54 

59 

60 

55 

1 


0 0 1 
0 0 1 
0 0 1 
0 0 1 


43 


ERROR MESSAGES 


The program gives two types of error messages. One is related to the 
problem size and the other is related to the finite-element (FE) modelling. 
In both cases the program STOPS normally and the error is flagged. The 
first type of error is easy to correct. Simply replace the dimension of the 
variable flagged with the new dimension value. This correction should be 
made in all of the parameter statements where that particular variable 
appears. These errors will pertain to the eight variables listed in the 
section Program Size Control. The second type of error is corrected by 
checking the FE model and the plastic analysis strategy. The three error 
messages and the probable corrections (remedies) are given below: 


1. ELEMENT ie HAD NONPOSITIVE DIAGONAL STIFFNESS MATRIX; ELEMENT 
CONNECTIVITY: i , j,k,l ,m,n,p,q. 

Check element "ie" connectivity which must follow the 
right-hand screw rule. 

2. IERR=I , NONPOSITIVE DEFINITE MATRIX 

Check boundary conditions to suppress all six rigid body 
modes of deformations. Impose restraints such that 
the three translations and the three rotations about 
x, y, and z axes are suppressed. 

3. NO CONVERGENCE AFTER 100 ITERATIONS 

This is the most difficult error to diagnosis and correct, but 
use the following guidelines. Increase the value of MAXIT, 
decrease the load-increment size, or increase the allowable error 
(ERIT) for stress convergence. Also, do not use a fine mesh 
(small elements) at a point load. If a fine mesh is used, then 
make the element immediately under the point load elastic. 
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IV. EXAMPLE PROBLEMS 


ELASTIC ANALYSIS OF SURFACE CRACK IN A PLATE 
A semi -ell iptical surface crack emanating from the root of a single- 
edge-notched plate subjected to remote stress (S = 1) was analyzed. The 
plate width was W = 45t, notch radius was r - 3t, and the thickness was t = 
1.0. The aspect ratio of the crack was a/c =■ 2.0, and the ratio of a/t = 
0.8. Because of symmetry, only one-quarter of the plate was modelled using 
8-node isoparametric elements. The model had 5380 nodes and 4375 elements. 
The FE idealization near the notch root is shown in the Figure 10(a) and the 
FE idealization at the crack plane is shown in the Figure 10(b). The crack 
front is divided into thirteen segments with larger segment lengths near the 
x.-axis. The mesh is nearly orthogonal at the deep section of the crack 
front and skewed towards the free surface (intersection of the crack front 
and the notch root). The domain used for the J-integral is shown in the 
Figure 10(c). The elastic modulus and Poisson's ratio were 1.0 and 0.3, 
respectively. In this analysis, N0DE1 was specified to be zero to eliminate 
printout of nodal displacements, forces and stresses. Fracture mechanics 
parameters G and J are calculated along the crack front from the midplane to 
the free surface. Although the mesh is not exactly orthogonal, both G 
calculated from the virtual -crack-closure technique (VCCT) and J calculated 
from the equivalent-domain integral (EDI) method agree well. The J-integral 
is not calculated where the crack front intersects the free surface (y-axis) 
for reasons explained in Reference 14. Three files, namely, surac2.dat 
(analysis file), SURAC2M (mesh file - only a partial listing), and the 
output file (surac2.res) are listed separately. All three files are 
provided along with the program file. 
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Crack front 
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(a) Finite-element idealization (c) Typical section across the crack front 

Figure 10. Finite-element idealization of a single edge notched tension specimen 



ELASTIC-PLASTIC ANALYSIS OF SINGLE-EDGE-CRACKED 
PLATE SUBJECTED TO TENSION 


A single-edge-crack-tension specimen (see Fig. 11(a)) with a = 0.5m, W 
= 1.0m, B = 0.5m, and H = 4.m, made of A36 steel with E = 207,000 MPa, v = 
0.3 and elastic limit stress a Q = 283.4 MPa was analyzed. The material 

stress- strain curve was represented by nine linear segments given in the 
analysis file "sent.dat". Figure 11(b) shows the FE model of one-fourth of 
the specimen. The model had five unequal layers in the z-direction, with 
total of 700 elements and 990 nodes. The nodal coordinates and the element 
connectivity are given in the file "SENTM" and a partial listing of the file 
is included here. The analysis data file is given in file "sent.dat" and a 
listing is also included in the manual. The first yield point occurred in 
element 693 at a stress level S = 19.15 MPa. The plastic analysis was 
continued to a value of S = 80 MPa. Three domain definitions (as shown in 
the Fig. 11(c)) were used to calculate the J-integral at several points 
along the crack front. In the three domains shown in Figure 11(c), the 
nodes on the outer boundary of the domain had deLorenzi's s-function equal 
to zero (s - 0) and all of the remaining nodes had s = 1. The J-integral s 
were calculated at 0.0 (midplane), 0.10, 0.18, 0.22, and 0.24mm from the 
midplane of the specimen. All fifteen values of the J-integral s, three at 
each point, are given in the output file. Listings of three files: the 
analysis file (sent.dat), the mesh file (SENTM), and the output file 
(sent. res) are given in the manual. All three files are included with the 
program file. 
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Figure 11.- Single-edge-crack- tension specimen and the finite-element model 


LISTING OF ANALYSIS, MESH AND OUTPUT FILES 
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INCREMENTAL LOAD FACTOR= 0.3000 
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SENTM 

MESH FILE: SENTM; 990-NODES, 700-ELEMENTS 
. OOOOOOOE+OO .2000000E+01 .OOOOOOOE+OO 

. OOOOOOOE+OO .2000000E+01 .lOOOOOOE+OO 

.OOOOOOOE+OO .2000000E+01 .1800000E+00 

.OOOOOOOE+OO .2000000E+01 .2200000E+00 

.OOOOOOOE+OO .2000000E+01 .2400000E+00 

.OOOOOOOE+OO .2000000E+01 .2500000E+00 

. 1200000E+00 .2000000E+01 .OOOOOOOE+OO 

. 1200000E+00 .2000000E+01 .1000000E+00 

. 1200000E+00 . 2000000E+01 .1800000E+00 

. 1200000E+00 . 2000000E+01 .2200000E+00 

. 1200000E+00 . 2000000E+01 .2400000E+00 

. 1200000E+00 . 2000000E+01 .2500000E+00 

. 2200000E+00 . 2000000E+01 .OOOOOOOE+OO 

.2200000E+00 . 2000000E+01 .lOOOOOOE+OO 

.2200000E+00 . 2000000E+01 .1800000E+00 

. 2200000E+00 . 2000000E+01 .2200000E+00 

.2200000E+00 . 2000000E+01 .2400000E+00 

.2200000E+00 . 2000000E+01 .2500000E+00 

. 3200000E+00 . 2000000E+01 .OOOOOOOE+OO 

. 3200000E+00 . 2000000E+01 .lOOOOOOE+OO 

. 3200000E+00 . 2000000E+01 .1800000E+00 

. 3200000E+00 . 2000000E+01 .2200000E+00 

. 3200000E+00 . 2000000E+01 .2400000E+00 

. 3200000E+00 . 2000000E+01 .2500000E+00 

. 3950000E+00 . 2000000E+01 .OOOOOOOE+OO 


*** CONTINUED *** 


966 

. 6050000E+00 

.OOOOOOOE+OO 

. 2500000E+00 

967 

. 6800000E+00 

.OOOOOOOE+OO 

.OOOOOOOE+OO 

968 

. 6800000E+00 

.OOOOOOOE+OO 

. lOOOOOOE+OO 

969 

. 6800000E+00 

.OOOOOOOE+OO 

. 1800000E+00 

970 

. 6800000E+00 

.OOOOOOOE+OO 

. 2200000E+00 

971 

. 680000 0E+00 

.OOOOOOOE+OO 

. 2400000E+00 

972 

. 6800000E+00 

.OOOOOOOE+OO 

. 2500000E+00 

973 

. 7800000E+00 

.OOOOOOOE+OO 

.OOOOOOOE+OO 

974 

. 7800000E+00 

.OOOOOOOE+OO 

. lOOOOOOE+OO 

975 

. 7800 000E+00 

.OOOOOOOE+OO 

. 1800000E+00 

976 

. 7800000E+00 

.OOOOOOOE+OO 

. 2200000E+00 

977 

. 7 8 00000E+00 

.OOOOOOOE+OO 

. 2400000E+00 

978 

. 7800000E+00 

.OOOOOOOE+OO 

. 2500000E+00 
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. 8 8 00000E+00 

.OOOOOOOE+OO 
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. 8800000E+00 
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. 1000000E+01 
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989 

. 1000000E+01 

.OOOOOOOE+OO 

. 2400000E+00 

990 

. 1000000E+01 

.OOOOOOOE+OO 

. 2500000E+00 
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26 
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SENTM 
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*** CONTINUED *** 
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1 
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No. OF SYMMETRIC BOUNDARY CONDITIONS 
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NUMBER OF G-CALCULATIONS = 5 

CRACK PLANE SYMMETRY Y/N 1 , 2 , 3 / 0 : 

G-REGI 0 N NO . = 1 

NO. OF ELEMENTS = 2 

ELEMENT NOS: 133 134 
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Plastic load vector time, Secs = 1.1879 Wall* 9.8317 

Total Time Used in CONTACT 0.8493E-04 

Solution Converged : 14 Iterations. Load = 0.59377E402 
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